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Computational  Analyses  and  Simulation  of  Toxic  Particle  Deposition 
in  the  Human  Respiratory  System 


1.  INTRODUCTION 

As  recently  reviewed  (see  COT/NA,  2003),  jet  fuel  performance  enhancers,  i.e.,  additives 
such  as  BHT,  DiEGME,  EGME,  BTEX,  etc.,  are  highly  toxic  and  may  pose  a  significant  health 
risk  to  Air  Force  personnel,  fuelers,  and  people  living  near  military  fields.  In  order  to  understand 
the  possible  health  effects  from  fuel  exposure  in  form  of  inhaled  droplets  and  vapor,  it  is  crucial 
to  know  where  and  at  what  concentrations  fuel  aerosols  and  their  vapors  deposit  in  the  human 
respiratory  system,  given  a  set  of  realistic  inlet  conditions. 

Direct  deposition  tests  of  toxic  fuel  in  human  lung  airways  are  impossible,  and  usually 
cost-intensive  experimental  deposition  data  lack  any  detailed  resolution.  Thus,  computational 
fluid-particle  dynamics  (CFPD)  simulations  offer  a  non-invasive,  less  expensive,  and  effective 
means  of  obtaining  fuel  vapor  and  droplet  transport  and  deposition  data  for  representative 
airways.  The  human  respiratory  system  has  very  complicated  structures.  The  oral  and  nasal 
airways  are  two  aerosol  entries  into  the  trachea.  After  the  trachea,  the  tracheobronchial  airways 
may  be  approximated  as  a  network  of  repeatedly  bifurcating  tubes  with  progressively  decreasing 
dimensions  (Weibel,  1963).  So  far,  computational  airflow  and  particle  transport  simulations  have 
concentrated  on  selected  parts  of  the  respiratory  system.  CFPD  studies  of  micro-size  particle 
deposition  in  human  oral  airways  and  tracheobronchial  airways,  employing  single,  double  and 
triple  bifurcating  models  have  been  extensively  conducted  by  several  research  groups  (Balashazy 
et  al.,  1999;  Comer  et  al.,  2001a,  b;  Kleinstreuer  and  Zhang,  2003;  Zhang  et  al.,  2002a,  b;  among 
others).  These  studies  assumed  that  the  size  of  inhaled  particles  stays  constant  in  the  respiratory 
tract.  However,  in  many  situations,  solid  or  liquid  particles  undergo  size  changes  because  they 
absorb  moisture  or  lose  mass  due  to  hygroscopity  or  evaporation.  The  effect  of  hygroscopity  on 
deposition  fraction  in  the  human  respiratory  system  has  been  investigated  experimentally  and 
theoretically  (Broday  and  Georgopoulos,  2001;  Ferron  et  al.,  1988;  1989;  Finlay  and  Stapleton, 
1995;  Gebhart  et  al.,  1990;  Hickey  and  Martonen,  1991;  Morrow,  1986;  among  others).  For 
example,  hygroscopity  can  substantially  increase  the  size  and  hence  the  airway  deposition 
fraction  of  drug  particles  with  initial  dry  sizes  of  about  0.5-2|i.m  (Ferron  et  al.,  1989).  In  contrast, 
as  highly  volatized  compounds,  jet  fuel  droplets  vaporize  quickly  and  may  undergo  dramatic  size 
reduction  as  they  move  through  the  respiratory  tract. 

Focusing  in  this  report  on  JP-8  fuel  droplet  and  vapor  deposition  in  a  human  upper  airway 
model,  the  contrasting  deposition  pattern  between  aerosols  and  vapor  as  well  as  the  resultant 
exposure  risk  are  of  interest.  Considering  different  inspiratory  flow  rates,  droplet  trajectories, 
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evaporation  and  deposition  are  modeled  with  Newton’s  Second  Law  and  scalar  equations,  i.e., 
using  a  decoupled  Euler-Lagrange  approach,  whereas  fuel  vapor  transport  and  deposition  is 
described  with  a  modified  mass  transfer  equation,  i.e.,  in  an  Eulerian  framework  (see 
Kleinstreuer,  2003). 

2.  THEORY 

2.1  Airway  Geometry 

The  detailed  descriptions  about  the  present  upper  airway  model  are  given  in  Zhang  & 
Kleinstreuer  (2003b),  which  includes  oral  cavity,  pharynx,  larynx  and  trachea,  and  a  symmetric 
triple  bifurcation  lung  airway  representing  generations  GO  (trachea)  to  G3  after  Weibel  (1963). 

2.2  Governing  Equations 

2.2.1  Airflow  and  heat  transfer 

In  order  to  capture  the  air  flow  structures  in  the  laminar-to-turbulent  flow  regimes,  i.e., 
300<Reiocai<104  for  the  present  airway  configuration  and  a  wide  range  of  inhalation  flow  rates 
(15<Qin  <60  l/min),  the  low-Reynolds-number  (LRN)  k-co  model  of  Wilcox  (1993)  was  selected 
and  adapted  which  has  been  demonstrated  to  be  appropriate  for  such  internal  flows  (Zhang  & 
Kleinstreuer,  2003a).  As  demonstrated  by  Zhang  et  al.  (2002a,  b),  the  transient  fluid-particle 
dynamics  can  be  accurately  matched  assuming  steady  flow,  i.e.,  using  equivalent  Reynolds-and- 
Stokes-number  pairs  for  all  inhalation  flow  rates.  All  air  flow  equations,  as  well  as  initial  and 
boundary  conditions,  are  given  in  Zhang  &  Kleinstreuer  (2003a,  b). 

2.2.2  Transport  of  micro-size  droplets 

With  any  given  ambient  concentration  of  non-interacting  spherical  droplets,  a  Lagrangian 
frame  of  reference  for  the  trajectory  computations  of  the  evaporating  droplets  can  be  employed. 
The  detailed  descriptions  of  micro-particle  trajectory  equation  with  turbulent  dispersion  are 
given  in  Zhang  et  al.  (2005). 

2.2.3  Droplet  evaporation 

This  analysis  considers  that  the  droplet  temperature  and  composition  are  maintained 
spatially  uniform  but  evolve  with  time  and  the  droplet  remains  spherical  during  vaporization. 
This  assumption  is  reasonable  because  of  the  small  Biot  number  (Bi«l)  for  the  micro-size 
droplets.  Considering  the  convective  heat  and  mass  transfer  over  the  whole  surface  of  the 
spherical  droplets,  the  change  in  droplet  mass  (or  size)  can  then  be  calculated  (see  Zhang  et  al., 
2004). 
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2.2.4  Mass  transfer  of  vapor  or  ultrafine  particles. 

The  convection  mass  transfer  equation  of  ultrafine  particles,  or  (JP-8)  fuel  vapor,  whose 
dominant  radial  transfer  mechanisms  are  Brownian  motion  and  turbulent  dispersion  can  be  given 
as  (Zhang  &  Kleinstreuer,  2003b,  2004): 


dY  o  /  v\ 


dt 


dxj 


dx, 


y  ' 

D  +  ^ 


V 


'Y  J 


dY_ 

dx 


(1) 


where  Y  is  the  mass  fraction,  <rY  =0.9  is  the  turbulence  Prandtl  number  for  Y.  JP-8  fuel  is  a 
multi-component  mixture  which  consists  mainly  of  C8  to  C16  paraffinic  hydrocarbons  with  other 
hydrocarbons  and  additives  also  present;  however,  the  diffusivity  in  air  does  not  vary 
significantly  from  compound  to  compound  (Gustafson  et  al.,  1997).  Thus,  the  conservative 
assumption  presently  implemented  was  to  set  D  =  0.05  cm2/sec  for  all  fractions  of  JP-8  jet  fuel  at 
T=293K  with  little  loss  in  accuracy  according  to  Gustafson  et  al.  (1997).  The  effect  of 
temperature  change  on  diffusivity  is  expressed  by  the  following  semi-empirical  correlation 
(Bejan,  1995): 

D(T)/D(T  =  293K)  =  (T /293)1'75  (2) 

The  aerosol  diffusion  coefficient  is  calculated  as  follows  (Finlay,  2001): 


Dp=(kBTCslip)/(  3>qidp)  (3) 

where  kB  is  the  Boltzmann  constant  (1.38xl0'23  J  K'1);  and  Cshp  is  the  Cunningham  slip 
correction  factor. 

Assuming  that  the  airway  wall  is  a  perfect  sink  for  aerosols  or  vapors  upon  touch,  the 
boundary  condition  on  the  wall  is  Yw=0.  This  assumption  is  reasonable  for  fast  gas-wall  reaction 
kinetics  (Fan  et  al.,  1996),  or  vapors  of  high  solubility  and  reactivity,  and  also  suitable  for 
estimating  the  maximum  deposition  of  toxic  vapor  in  the  airways.  For  less  soluble  vapors,  the 
wall  concentration  would  be  greater  than  zero  so  that  transport  in  tissue  and  in  airways  must  be 
considered  simultaneously  for  simulating  vapor  uptake.  Assuming  that  the  surface  of  respiratory 
epithelium  is  covered  by  a  mucus  layer  with  uniform  thickness  and  a  lipid  layer  lie  below  the 
mucus  lining  simulating  the  transport  barrier  by  the  epithelial  cell  membrane,  Keyhani  et  al. 
(1997)  derived  a  flux  condition  at  the  airway  boundary  including  the  vapor  transport  in  the 
tissue.  This  boundary  condition  was  inferred  from  the  mass  conservation  and  mass  diffusion  of 
vapor  from  the  airway  to  the  mucus  layer  and  to  the  tissue,  which  is  given  as: 

Din?f  +  KY  =  0  (4) 

on 

where  Din  is  the  inlet  tube  diameter,  n  is  the  direction  normal  to  the  airway  wall,  and  K  is  a 
dimensionless  parameter,  here  we  call  it  absorption  parameter,  which  is  defined  as 
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where  Da,Dm  are  the  vapor  diffusivity  in  the  air  and  the  liquid  mucus  phase,  respectively;  Hm  is 
the  thickness  of  the  mucus  layer;  and  ft  is  the  equilibrium  partition  coefficient  for  a  given 
contaminant  molecular,  which  can  be  determined  by  Henry”s  law;  %  =  fkr  /  Dm  with  kr  being  a 
single  rate  constant  considering  the  chemical  reactions  of  vapors  in  the  mucus  layer; 
Pt  =  ')DI  /  H,  is  the  lipid  permeability  coefficient  with  y,  D, ,  Hl  being  lipid-mucus  partition 


coefficient,  vapor  diffusivity  in  the  lipid  and  the  thickness  of  the  lipid  layer,  respectively. 

Clearly,  for  a  highly  soluble  (/?  -»  0)or  high  reactive  (  kr  — >  ■»)  vapor,  the  boundary 
condition  (4)  reduces  to  Y=0.  If  vapor  is  insoluble  (/?—>=*=),  the  boundary  condition  reduces  to 
the  zero  mass  flux  condition,  i.e.,  dY  ldn  =  0 .  Ignoring  removal  of  contaminant  molecules  in 
the  mucus  by  chemical  reaction  and  the  resistance  to  transport  across  the  lipid  barrier,  the 
absorption  parameter  K  reduces  to: 
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2.3  Deposition  Parameters 

The  regional  deposition  of  micro-droplets  in  human  airways  can  be  quantified  in  terms  of 
the  deposition  fraction  (DF)  or  deposition  efficiency  (DE)  in  a  specific  region  (e.g.  oral  airway, 
first,  second  and  third  bifurcations  etc.);  they  are  defined  as: 


DF, 


particle 


Number  of  deposited  particles  in  a  specific  region 
Number  of  particles  entering  the  mouth 
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DE 


particle 


Number  of  deposited  particles  in  a  specifc  region 
Number  of  particles  entering  this  region 


(8) 


The  DEs  and  DFs  are  the  same  for  the  oral  airway  model  in  this  study. 

The  deposition  fraction  (DF)  of  vapors  or  ultrafine  particles  can  be  calculated  with  the 
regional  mass  balance  or  the  sum  of  local  wall  mass  flux.  According  to  the  mass  balance  in  one 
specific  region,  DF  is  defined  as: 


where  A  is  the  tube  cross  sectional  area,  and  u  is  the  axial  velocity.  The  subscripts  “0”  and  “i” 
denote  the  properties  at  the  inlet  and  the  ith  daughter  tube  of  a  selected  generation,  respectively. 
In  the  oral  airway  model,  the  “i  ”  denotes  the  trachea  and  n=l. 

As  for  the  local  wall  mass  flux  of  vapors,  it  can  be  determined  as 

mw  =  -p4(£>  +  — )^|,  (10) 

aY  on 

where  A*  is  the  area  of  local  wall  cell  (i).  The  local  vapor  deposition  fraction,  which  is  defined 
as  the  ratio  of  local  wall  mass  flux  to  the  inlet  mass  flux,  can  be  expressed  as 


DFM  =[-4(S+^)f£|,]/(&,n.)  (11) 

Gy  on 

and  the  regional  deposition  fraction  can  be  determined  as 

n  y  ^  Y 

sf„.  =  lK(o*7-)yi«A)  (12) 

where  n  is  the  number  of  wall  cells  in  one  specific  airway  region,  e.g.,  oral  airway,  first  airway 
bifurcation,  etc.  The  local  vapor  deposition  patterns  can  be  quantified  in  terms  of  a  deposition 
enhancement  factor  (DEF)  (cf.  Balashazy  et  al.,  1999,  2003),  which  is  defined  as  the  ratio  of 
local  to  average  deposition  densities,  i.e., 


dn 

DEF  indicates  vapor  deposition  “hot  spots”  in  a  given  region. 
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2.4  Numerical  Method 

The  numerical  solutions  of  the  governing  equations  including  flow,  heat  and  mass 
transfer  equations  were  carried  out  with  a  user-enhanced  finite-volume  based  program,  i.e., 
CFX4.4  from  ANSYS,  Inc.,  previously  AEA  Technology  (2001).  The  numerical  program  uses  a 
structured,  multiblock,  body-fitted  coordinate  discretization  scheme  (Zhang  &  Kleinstreuer, 
2002).  A  Higher-Order  Upwind  (HUW)  differencing  scheme,  which  is  second-order  accurate  in 
space,  was  used  to  model  the  advective  terms  of  the  transport  equations.  The  sets  of  linearized 
and  discretized  equations  for  all  variables  were  solved  using  the  Block  Stone’s  method. 

The  droplet  transport  and  evaporation  equations  were  also  solved  with  the  commercial 
software  in  combination  with  user-supplied  programs.  The  particle  trajectory  equations  were 
solved  using  Gear’s  BDF  method.  After  each  iteration  for  each  particle,  the  information  about 
position,  time,  diameter,  temperature  and  three  components  of  the  velocity  as  well  as  the  speed 
with  which  the  particles  cross  the  control  volume  boundaries  was  obtained.  Droplet  deposition 
occurs  when  its  center  comes  within  a  radius  from  the  wall,  i.e.,  local  surface  effects  such  as 
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droplet  evaporation  (blowing),  migration,  or  resuspension  have  been  currently  ignored.  In  the 
present  model  the  number  of  particles,  n=l  0,000,  was  determined  by  increasing  the  inlet  particle 
concentration  until  the  deposition  fraction  became  independent  of  the  number  of  particles 
simulated. 

The  computational  mesh  was  generated  with  CFX  Build4.  The  near-wall  region  required 
a  very  dense  mesh.  Specifically,  the  thickness  of  the  near-wall  cells  was  chosen  to  fully  contain 
the  viscous  sub-layers  and  to  resolve  any  geometric  features  present  there.  The  mesh  topology 
was  determined  by  refining  the  mesh  until  grid  independence  of  the  flow  field  solution  and 
particle  deposition  fractions  was  achieved.  The  final  mesh  contains  about  420,000  and  670,000 
cells  for  the  oral  airway  and  four-generation  airway  model,  respectively.  The  computations  were 
performed  on  an  SGI  Origin  2400  workstation  with  32GB  RAM  and  multiple  450  MHz  CPUs. 
The  steady-state  solution  of  the  flow  field  was  assumed  to  be  converged  when  the  dimensionless 
mass  residual,  (Total  Mass  Residual)/(Mass  Flow  Rate)  <  10'3.  The  convergence  of  k  and  co  was 
monitored  as  well.  Typical  run  time  for  the  fluid  flow,  energy  and  mass  transfer  simulations  on 
eight  processors  with  parallel  algorithm  was  approximately  24  hours  for  the  oral  airway  model 
and  8  hours  for  the  four-generation  model.  Utilizing  the  converged  flow  field  solution,  the 
droplet  trajectory  simulations  required  approximately  4  to  6  hours  for  each  case  considered. 

3.  RESULTS 

3.1  Model  Validations 

For  accurate  computational  fluid-particle  dynamics  (CFPD)  simulations,  matching 
comparisons  with  theoretical  or  experimental  fluid  velocity  profiles,  pressure  drops  and  local 
particle  distributions  are  necessary.  Our  CFPD  model  has  been  validated  with  various 
experimental  data  sets  for  steady  and  transient  laminar  flows  in  bifurcations  (Comer  et  al,  2001a; 
Zhang  &  Kleinstreuer,  2002)  and  for  laminar,  transitional  and  turbulent  flows  in  tubes  with  local 
obstructions  (Kleinstreuer  &  Zhang,  2003;  Zhang  &  Kleinstreuer,  2003a).  Especially,  the  low- 
Reynolds-number  (LRN)  k-omega  model  has  been  extensively  validated  and  has  been  proven  to 
be  an  applicable  approach  to  capture  the  velocity  profiles  and  turbulence  kinetic  energy  for 
laminar-transitional-turbulent  flows  in  the  constricted  tubes  of  the  upper  airways  (see  Zhang  & 
Kleinstreuer,  2003a).  Similarly,  the  simulated  spherical,  micro-particle  depositions  in  airways 
were  successfully  compared  with  measured  deposition  efficiencies  and  deposition  patterns 
(Comer  et  al.,  2001b;  Zhang  et  al.,  2002a-d).  The  simulation  of  spherical,  nanoparticle  deposition 
due  to  diffusional  transport  has  been  validated  with  both  analytical  solutions  in  straight  pipes  and 
experimental  data  for  a  double-bifurcation  airway  model  (Shi  et  al.,  2004)  as  well  as 
experimental  data  in  an  oral  airway  model  (Zhang  &  Kleinstreuer,  2003b).  For  example,  particle 
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deposition  fractions  as  a  function  of  Stokes  number  are  compared  with  the  observations  by 
Cheng  et  al.  (1999)  in  Fig.  1  for  three  inhalation  rates.  The  Stokes  number  is  defined  here  as  by 

Cheng  et  al.  (1999),  i.e.,  St=ppdp2  U/9pD,  with  pp  being  the  particle  density,  dp  being  the 

particle  diameter,  and  U  being  the  mean  velocity  evaluated  as  (Q/A),  where  A  is  the  mean  cross- 
sectional  area,  and  D  is  the  minimum  hydraulic  diameter.  The  validation  result  for  droplet 
vaporization  simulations  is  given  in  Fig.  2.  The  simulated  temporal  diameter  variation  of  a  JP-8 
fuel  droplet  is  compared  with  experiment  measurements  conducted  by  Runge  et  al.  (1998), 
where  a  JP-8  fuel  droplet  with  an  initial  diameter  of  639pm  was  suspended  in  294K  air  moving 
at  a  constant  speed  of  3m/s.  It  can  be  seen  that  the  simulated  vaporization  law  for  JP-8  fuel,  using 
A#e£7;Ve=0.05cm2/sec,  is  in  excellent  agreement  with  experimental  data. 


Stokes  Number,  St=ppdpU/(9jiD) 

Fig.  1:  Comparison  of  simulated  fine  particle  deposition  fractions  in  the  oral  airway  model  with  the 
experimental  data  of  Cheng  et  al.  (1999). 
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Fig.  2:  Comparison  of  diameter  evolution  for  a  JP-8  droplet  with  the  experimental  data 

of  Runge  et  al  (1998). 
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In  summary,  the  good  agreements  between  experimental  findings  and  theoretical 
predictions  instill  confidence  that  the  present  computer  simulation  model  is  sufficiently  accurate 
to  analyze  transport  and  deposition  of  toxic  particle/droplet  and  vapor  in  three-dimensional  oral 
and  upper  bronchial  airways  associated  with  laminar-to-turbulent  airflows. 

3.2  Airflow  Fields 

Some  selected  results  from  our  simulations  for  different  inspiratory  flow  rates  and  inlet 
air  temperatures  are  given  below.  Figures  3a  and  3b  show  mean  velocity  profiles  in  an  oral 
airway  model  starting  from  mouth  to  trachea  and  a  bifurcating  airway  model  GO  to  G3  (GO  is  the 
trachea)  with  an  inspiratory  flow  rate  of  Qin=30  1/min  and  Tjn=310K,  respectively.  The  selected 
cross-sectional  views  display  the  axial  velocity  contours  as  well  as  secondary  velocity  vectors. 
Clearly,  skewed  velocity  profiles  generated  by  the  centrifugal  force  can  be  observed  in  the 
curved  portion  from  the  oral  cavity  to  the  pharynx/larynx.  Moreover,  a  central  asymmetric  jet 
and  a  recirculation  zone  are  created  in  the  laryngeal  region  because  of  the  restriction  of  the 
(assumed  stationary)  vocal  folds.  The  secondary  motion  is  set  up  when  the  flow  turns  a  bend 
from  the  mouth  to  the  pharynx  because  of  the  centrifugally  induced  pressure  gradient  (see 
Kleinstreuer  &  Zhang,  2003).  Turbulence  may  locally  occur  during  this  medium-level  breathing 
(Qin=30  1/min),  where  the  turbulent  fluctuations  are  weak  in  the  oral  cavity;  but,  they  becomes 
strong  after  the  constriction  of  soft  palate  and  rises  rapidly  after  the  glottis,  eventually  decaying 
more  slowly  while  approaching  an  asymptotic  level.  With  further  redistribution  of  the  kinetic 
energy  of  the  flow  over  most  of  the  cross  section  accompanied  by  the  onset  of  turbulence,  the 
velocity  profiles  become  more  blunt  from  the  left  to  the  right  wall,  and  the  maximum  velocity 
zone  moves  to  the  anterior  wall  at  the  cross  section  six  diameters  from  the  glottis  (see  Fig.  3a). 
Fig.  3b  shows  the  mean-velocity  fields  in  the  planar  triple-bifurcation  airways  which  are 
physiologically  90°  turned  with  respect  to  the  oral  airways.  The  velocity  profiles  and  turbulence 
quantities  extracted  from  the  trachea  in  the  oral  airway  model  were  adjusted  as  the  inlet 
conditions  for  the  bifurcating  airway  model  GO  to  G3.  The  air  stream  in  the  branching  airways 
splits  at  each  flow  divider  and  new  boundary  layers  are  generated  at  the  inner  walls  of  the 
daughter  tubes.  The  velocity  profiles  are  naturally  skewed  in  each  daughter  tube,  and  hence,  each 
daughter  tube,  or  generation,  may  experience  a  different  flow  rate.  The  primary  characteristics  of 
the  flow  fields  at  cross  sections  C-C',  D-D',  and  E-E'  are  quite  similar;  specifically:  (i)  skewed 
velocity  profiles  with  the  maximum  velocity  near  the  inner  wall  around  the  divider;  and  (ii)  two 
distinct  secondary  vortices  appearing  at  the  upper  and  lower  side  of  the  tube,  which  moves  the 
high  speed  fluid  up  around  the  top  of  the  tube  towards  the  outside  of  the  bifurcation  and  low 
speed  fluid  from  the  outside  of  the  bifurcation  along  the  symmetry  plane  towards  the  inside  of 
the  bifurcation. 
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Fig-3  :  Side-plane  or  mid-plane  velocity  profiles  as  well  as  selected  cross-sectional  axial  velocity 
contours  and  secondary  velocity  vectors  with  an  inspiratory  flow  rate  Qin=30  1/min  in:  (a)  an  oral  airway 
model;  and  (b)  a  bifurcating  airway  model  GO  to  G3. 


The  details  about  the  airflow  structures  and  particle  distributions  in  both  oral  and 
tracheobronchial  airways  under  cyclic  breathing  conditions  as  well  as  the  effects  of  airway 
geometric  features,  inspiratory  flow  rate  and  thermal  condition  of  incoming  air  can  be  found  in 
Zhang  &  Kleinstreuer  (2002,  2003b),  Zhang  et  al.  (2002d,  2004). 


3.3  JP-8  Fuel  Aerosols 


Inertial  impaction  is  the  main  deposition  mechanism  for  micron-size  droplets;  hence, 
particle  size  greatly  affects  droplet  deposition.  However,  jet  fuel  is  a  multi-component  liquid  and 
the  C8  to  Cll  fractions  are  especially  volatile.  Here,  one  extreme  case  is  presented,  i.e.,  it  is 
assumed  that  the  inlet  vapor  fraction  of  JP-8  fuel  is  zero  and  the  background  concentration  of 
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vapor,  created  by  droplet  evaporation,  is  negligible.  Clearly,  the  diameter  of  each  droplet 
decreases  gradually  in  the  airways  due  to  the  loss  of  its  mass  by  evaporation  which  is  a  function 
of  the  airflow  momentum  as  well  as  the  droplet  heat  and  mass  transfer.  Affected  by  the  highly 
non-uniform,  occasionally  turbulent  airflow  structures  in  the  oral  airways,  the  trajectories  for 
different  droplets  vary;  hence,  the  changes  of  droplet  diameters  are  different  because  of  the 
different  local  airflow  as  well  as  heat  and  mass  transport.  Thus,  droplets  enter  the  mouth 
monodisperse  but  turn  polydisperse  during  vaporization.  Figure  4  shows  the  cumulative 
distribution  function  of  the  droplet  size  at  the  outlet  of  the  oral  airways,  or  the  inlet  to  the  first 
bifurcation,  for  different  inhalation  rates  and  inlet  air  temperatures.  The  fraction  less  than  a 
particular  size  can  be  obtained  directly  from  the  graph.  The  fraction  of  particles  having  diameters 
between  two  sizes  can  be  determined  by  the  difference  of  the  cumulative  fractions  at  these  two 

sizes.  The  initial  droplet  size  is  a  uniform  distribution  with  dp  =  5 /An .  As  expected,  the  change 

of  droplet  size  decreases  with  increasing  flow  rate  due  to  decreasing  residence  times.  The 
variation  of  the  droplet  size  is  larger  with  higher  inlet  air  temperature  because  the  saturation 
vapor  pressure  at  the  liquid/gas  interface  increases  with  temperature.  The  mean  droplet  diameters 
after  the  oral  airways  are  dp=4.03,  4.57,  4.83  pm  for  Qin=15,  30,  60  1/min  with  Tm=263K, 
respectively,  and  dp=2.80  pm  for  Qin=  30  1/min  with  Tin=310K.  The  variations  in  droplet  size 
should  greatly  influence  droplet  deposition. 


Fig.  4:  Cumulative  distribution  of  JP-8  fuel  droplet  sizes  at  the  inlet 
of  the  bifurcation  airway  model. 
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The  local  deposition  patterns  in  the  oral  airway  model  are  shown  in  Fig.  5  for  initially 
dp=5|im  and  Qi„=30  1/min  with  and  without  evaporation.  Due  to  inertial  impaction,  droplets 
mainly  deposit  at  stagnation  points  for  axial  particle  motion,  such  as  the  tongue  portion  in  the 
oral  cavity,  the  outer  bend  of  the  pharynx/larynx,  and  the  regions  just  upstream  of  the  glottis.  A 
few  droplets  may  also  deposit  outside  these  pronounced  regions,  being  influenced  by  turbulent 
dispersion  as  well  as  secondary  and  recirculating  flows.  Clearly,  droplet  deposition  decreases  due 
to  the  reduced  droplet  diameters  with  evaporation. 


Qln=30  l/min,  dp=5pm 
(with  evaporation) 


Qln=30  l/mln,  dp=5p.m 
(without  evaporation) 


Fig.  5:  Droplet  deposition  patterns  in  the  oral  airway  model. 

Droplet  evaporation  is  very  important  in  simulating  the  deposition  of  volatile  aerosols, 
e.g.,  JP-8  fuel  droplets,  in  human  airways.  This  also  can  be  seen  from  the  comparisons  of 
predicted  droplet  deposition  fractions  (DF)  in  the  oral  airways  for  different  inhalation  conditions 
(see  Fig.  6).  Figure  6  indicates  that  the  DFs  of  droplets  decrease  significantly  if  vaporization 
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occurs.  The  DF  increases  with  increasing  initial  particle  diameter  and  inhalation  flow  rate  due  to 
enhanced  inertial  impaction  (Fig.  6). 


Droplet  diameter  (|xm) 

Fig.  6:  Deposition  fractions  of  JP-8  fuel  droplets  in  the  oral  airway  model. 

Figure  7  depicts  the  overall  droplet  deposition  efficiency  (DE)  in  the  upper  bronchial 
airways  (G0-G3),  which  is  defined  as  the  number  of  deposited  droplets  to  that  entering  at  the 
inlet  of  airway  segment  GO  to  G3.  The  deposition  efficiency  with  evaporation  (Tjn=263K)  is  still 
lower  than  that  without  evaporation.  There  is  no  droplet  deposition  for  the  case  of  dp=  3  pm  with 
evaporation.  The  percent  decrease  in  the  deposition  efficiency  due  to  evaporation  is  higher  for 
the  relatively  smaller  size  droplets  (say,  3  and  5pm).  This  is  because  large  size  particles  tend  to 
shrink  much  slower  than  smaller  particles.  The  overall  DE  in  GO  to  G3  with  an  ideal  inlet 
condition,  i.e.,  fully-developed  air  flow  with  corresponding  (parabolic)  particle  distribution,  is 
also  given  in  Fig.  7  for  comparison.  Clearly,  the  inlet  air  velocity  profile  and  particle  distribution 
strongly  affect  the  particle  deposition  efficiency,  where  idealized  inlet  flow  conditions  may 
overpredict  actual  deposition  efficiencies.  However,  regardless  of  the  simulated  cases,  micro-size 
particles  tend  to  deposit  at  the  carinal  ridges  in  the  bronchial  airways  due  to  inertia  impaction 
(see  Fig.  8  for  droplets  in  the  present  case,  and  Fig.  4  in  Zhang  et  al.  (2002b)  for  solid  particles). 
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Droplet  Size  (|itn) 

Fig.  7:  Deposition  efficiencies  of  JP-8  fuel  droplets  in  the  bifurcation  airway  model  with  an  inspiratory 
flow  rate  Qin=30  1/min.  The  inlet  air  temperature  is  263K  for  the  evaporation  case. 


Fig.  8:  Droplet  deposition  patterns  in  the  bifurcation  airway  model  (Qin=30  1/min). 

Droplet  diameter  is  7  pm  at  mouth  inlet  and  evaporation  is  considered. 

3.4  JP-8  Fuel  Vapor 

As  aforementioned,  the  degree  of  vapor  absorption  by  the  airway  wall  may  influence  fuel 
vapor  transport  and  deposition  in  the  airways.  Because  JP-8  fuel  is  a  multi-component  mixture 
and  the  solubility  of  its  components  changes  very  largely  (Gustafson  et  al.,  1997),  the  absorption 
parameter  K  in  Eq.  (5)  may  vary  from  10'3  to  103  for  different  compounds.  As  a  case  for 
estimating  the  maximum  deposition  of  toxic  vapor  in  the  airways,  the  following  results  are  based 
on  the  perfect  absorption  assumption  of  vapor  in  the  airway  surface  (i.e.,  YWaii=0),  except  in 
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Section  3.4.4  which  discusses  the  impacts  of  the  absorption  parameter  K. 


3.4.1  Distributions  of  mass  fractions  and  deposition  patterns 

Figures  9a  and  9b  display  the  concentration  contours  of  JP-8  fuel  vapors  in  the  oral 
airway  and  airway  segment  GO  -  G3.  A  uniform  inlet  JP-8  vapor  concentration  condition  was 


Fig.  9:  Concentration  contours  of  JP-8  vapor  in  at  Qjn=15  1/min  and  T;n=263K  in:  (a)  the  oral  airway 
model  (y=0  plane);  and  (b)  the  bifurcation  airway  model  (z=0  plane). 


assumed.  The  distributions  of  vapor  concentration  are  dependent  of  velocity  profiles  so  that 
areas  of  low  velocities  are  associated  with  low  concentrations  due  to  the  low  convective  mass 
transfer.  The  local  concentration  gradient  near  the  wall  indicates  the  vapor  deposition  rate  at  that 
site.  Thus,  the  local  deposition  rate  tends  to  be  high  at  the  oral  cavity  and  becomes  low  in  the 
trachea  due  to  the  vapor  redistribution  over  the  cross  sections  driven  by  secondary  and  reverse 
flows.  The  vapor  field  in  the  trachea  for  Qin=15  1/min  cannot  be  mixed  as  rapidly  as  the  heat, 
where  the  temperature  distributions  tend  to  be  uniform  at  E-E'  and  F-F'.  This  may  be  attributed 
to  the  relatively  high  Lewis  number  for  JP-8  vapor  (Le=4).  The  thermal  effect  on  concentration 
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distribution  is  apparent  for  this  low  inspiratory  flow  rate  case,  but  it  tends  to  be  negligible  for 
high  flow  rate  cases  (Qin=30  and  60  1/min)  (Zhang  &  Kleinstreuer,  2003b).  In  airway  segment 
G0-G3,  the  mass  transfer  patterns  are  also  strongly  affected  by  the  skewed,  asymmetric  velocity 
profiles  which  indicate  different  vapor  deposition  at  different  airway  tubes. 

Figure  10  shows  examples  of  vapor  deposition  patterns  in  terms  of  distributions  of  DEF 
(see  Eq.  (13))  in  the  oral  airway  model  and  the  bifurcation  airway  model.  Clearly,  the  deposition 
patterns  are  somewhat  inhomogeneous.  In  the  oral  airway  model,  the  enhanced  deposition  may 
occur  at  the  entrance,  outside  bend  of  the  pharynx  and  the  throat  because  of  great  degrees  of 
mixing,  i.e.,  large  concentration  gradients  at  these  regions.  Turning  to  the  bifurcation  airway 
model  GO  to  G3,  the  enhanced  deposition  mainly  occurs  at  the  cranial  ridges  and  the  inside  walls 
around  the  carinal  ridges  due  to  the  complicated  air  flows  and  large  particle  concentration 
gradients  in  these  regions.  Of  interest  is  that  the  maximum  DEF-values  for  vapor  are  much  lower 
than  those  for  micro-droplets,  i.e.,  the  DEFmax  for  micro-particles  is  of  the  order  of  10  to  10 
(see  Zhang  et  al.,  2005),  while  DEFmaX  for  vapor  is  of  the  order  of  1  (see  Fig.  10).  This  indicates 
that  the  deposited  vapor  is  much  more  uniformly  distributed  in  the  upper  airways  when 
compared  to  the  micro-droplet  deposition. 

(a)  (b) 


Fig.  10:  3-D  distributions  of  deposition  enhancement  factor  (DEF)  of  JP-8  fuel  vapor  under  steady 
inhalation  with  Qin=30  1/min  in:  (a)  the  oral  airway  model;  and  (b)  the  bifurcation  airway  model. 
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3.4.2  Deposition  fractions  of  fuel  vapor 

The  deposition  fractions  of  JP-8  fuel  vapor  in  the  upper  airways  under  different  inhalation 
conditions  are  shown  in  Fig.  lla-c.  It  can  be  seen  from  Fig.  11a  that  higher  vapor  deposition 
corresponds  to  larger  inlet  air  temperatures,  especially  for  the  high  inhalation  flow  rate  case.  This 
can  be  attributed  to  the  relatively  large  diffusion  coefficient  at  relatively  high  temperatures  (see 
Eq.  (2)).  Generally  speaking,  the  deposition  fraction  is  weakly  affected  by  the  inlet  air 
temperature  for  both  low  and  high  inspiratory  flow  rates.  In  other  words,  the  thermal  effects 
could  be  negligible  when  calculating  the  total  or  segmental  deposition  fraction  of  vapors  or 
ultrafme  particles  in  the  upper  airways.  The  flow  rate  also  shows  a  significant  effect  on  the 
deposition  of  vapor  (cf.  Fig.  1  lb  and  c).  The  higher  the  flow  rate,  the  lower  is  the  deposition 
fraction.  This  may  be  because  of  the  longer  residence  times  for  vapors  with  low  flow  rates  and  is 
consistent  with  the  experimental  observations  of  Li  et  al.  (1998)  for  deposition  of  Iodine  vapor  in 
a  tracheobronchial  cast.  The  deposition  fraction  may  increase  75%  to  120%  in  the  oral  airways 
and  GO  to  G3  when  switching  from  exercise  breathing  ( Qin  =60l/min)  to  low-level  breathing 
(Qjn  =  151/min).  The  deposition  fraction  of  JP-8  fuel  vapor  in  the  oral  airway  and  GO  to  G3  under 
resting  condition  can  be  as  high  as  36%  and  25%,  respectively. 


3.4.3  Calculation  of  mass  transfer  coefficient 

Calculation  of  the  respiratory  mass  transfer  coefficient,  hm ,  is  helpful  in  quantitatively 

predicting  the  regional  uptake  of  inhaled  vapor  (Cheng  et  al.,  1997a,  b).  The  mass  balance  for 
one  airway  unit  is: 

N-Amll  ^m,Y, -m0Y0  (14) 

1=1 


which  yields  the  wall  mass  flux  of  species  N ,  so  that 


pYm  ~pYwa„ 

For  simplicity,  we  only  used  the  average  values  of  the  inlet  and  outlet  cross-sectional  mass 
fractions  to  determine  the  regional  pYm  ,  i.e., 


(15) 


(16) 


The  Sherwood  number  (Sh),  characteristic  of  convective  dispersion,  can  be  calculated  as: 

Sh  =  (hmD)/D  (17) 

The  diameters  in  the  trachea  and  the  parent  tubes  are  employed  to  determine  the  Sh-number  for 
the  oral  airways  and  bifurcation  units,  respectively. 
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(a) 


Distance  from  the  mouth  (cm) 


(b) 


Fig.  11:  Deposition  fractions  of  JP-8  fuel  vapor  in  the  upper  human  airway 

Figure  12a  depicts  the  regional  Sherwood  number  (Sh)  for  the  oral  airway  model  vs.  the 
product  of  Reynolds  (Re)  and  Schmidt  (Sc)  numbers,  both  evaluated  in  the  trachea,  based  on 
different  inhalation  flow  rates  and  diffusion  coefficients  (i.e.,  particle  sizes).  The  correlation 
between  Sh-number  and  ReSc-group  can  be  expressed  as: 
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Sh  =  0.852(Re1S'c)0405  1500  <  Re  <  7000,  1  <  Sc  <  230  (18) 

where  r2=0.99.  The  convective  mass  transfer  coefficients  for  JP-8  fuel  vapor  in  terms  of  the  Sh- 
number  has  been  also  obtained  for  each  individual  bifurcation  in  airway  generations  GO  to  G3 
(Fig.  12b).  Affected  by  the  non-homogeneous  airflow  structures  and  mass  concentration 
distributions,  the  regional  Sh-numbers  are  slightly  different  for  each  bifurcation.  However,  a 
best-fit  correlation  yielded  (Fig.  12b): 

Sh  =  1.727(R e  Sc)0  328  600  <  Re  <  6000,  Sc  -  3  (19) 

where  ?=Q.ll. 


(a) 


Fig.  12:  Regional  Sherwood  number  (Sh)  versus  the  product  of  Re  and  Sc  for:  (a)  human  oral  airways; 
and  (b)  human  upper  bronchial  tree. 
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3.4.4  Impact  of  airway  wall  absorption  parameter  K 

The  absorption  parameter  K  for  different  compounds  of  JP-8  fuel  may  vary  from  10'3  to 
103,  resulting  in  large  variations  of  the  deposition  fractions.  Figure  13a  shows  the  variation  of 
deposition  fraction  in  the  oral  and  bifurcation  airway  models  as  a  function  of  K.  When  K  is  less 
than  one,  the  deposition  is  very  low  in  the  upper  airways  due  to  the  low  solubility  of  species  in 
the  mucus  layer.  The  deposition  fraction  is  greatly  dependent  on  K  for  1<K<1000.  If  K>1000, 
the  deposition  fraction  is  very  close  to  that  for  the  perfectly  absorbing  wall  condition  (i.e., 
Ywaii=0).  The  impact  of  K  on  the  species  mass  transfer  coefficients  in  the  airways  are  depicted  in 
Fig.  13b.  Clearly,  K  almost  has  no  effects  on  the  mass  transfer  coefficients  of  species  with  the 
same  diffusivities. 
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Fig.  13:  The  effects  of  airway  wall  absorption  parameter  K  on:  (a)  deposition  fraction  and  (b)  mass 
transfer  coefficients  in  the  human  upper  airway  model  with  Qin=30  1/min  and  Da  =0.05cm/s2. 
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4.  CONCLUSIONS 


4.1  Fundamentals  and  New  Discoveries 

The  first  part  of  the  three-year  project  focused  on  the  fundamentals  of  “lung  aerosol 
dynamics”  and  generated  the  following  results: 

(1)  Matching  Stokes  and  Reynolds  numbers  of  cyclic  inspiratory  waveforms,  i.e.,  the  average 
between  mean  and  peak  values,  was  found  to  generate  about  the  same  particle  deposition  under 
cyclic  flow  condition  as  steady-state  simulation  results.  Hence,  with  the  corresponding  matching 
Reynolds  and  Stokes  numbers,  the  extrapolation  of  constant-flow  deposition  data  to  cyclic-flow 
situation  at  low  Womersley  numbers  is  possible  (see  Zhang  et  al.,  2002a,  b). 

(2)  The  Low-Reynolds-number  (LRN)  k-(0  model  for  internal  flow  was  found  to  be  appropriate 
for  simulating  all  airflow  regimes,  i.e.,  laminar,  transitional  and  turbulent,  in  the  human  upper 
airways  (see  Kleinstreuer  &  Zhang,  2003;  Zhang  &  Kleinstreuer,  2003a;  Zhang  et  al.,  2002c). 

(3)  Depositions  of  both  micro-  and  nano-size  particles  vary  measurably  in  the  human  upper 
airways;  however,  the  deposition  distributions  are  much  more  uniform  for  nanoparticles.  The 
maximum  deposition  enhancement  factor,  which  is  defined  as  the  ratio  of  local  to  average 
deposition  concentrations,  ranges  from  about  40  to  2400  for  microparticles  and  about  2  to  1 1  for 
nanoparticles  with  inspiratory  flow  rates  in  the  range  of  15<Qjn<60  1/min.  The  nearly  uniform 
distribution  of  deposited  nanoparticles  may  imply  greater  toxicity  effects  of  such  particles  when 
compared  to  microparticles  made  of  the  same  material  (see  Zhang  &  Kleinstreuer,  2004;  Zhang 
et  al.,  2005). 

4.2  JP-8  Jet  Fuel  Droplets  and  Vapor 

The  transport  and  deposition  of  JP-8  jet  fuel  droplets  and  vapor  in  a  human  upper  airway 
model  were  simulated  and  analyzed  for  different  quasi-steady  inhalation  conditions.  The 
experimentally  validated  CFPD  results  show  the  following: 

(1)  Thermal  effects  in  the  oral  airways  are  considerable  at  low-level  breathing  (Qjn=  1 5  1/min), 
where  changing  temperature  distributions  influence  the  velocity  fields  measurably,  especially  in 
the  trachea.  However,  thermal  effects  tend  to  be  minor  for  medium  to  high-level  breathing 
(Qin=30  to  60  1/min).  Thermal  effects  in  the  bronchial  airways  are  minor  for  all  inhalation  flow 
rates  due  to  the  warming-up  of  the  air  when  moving  through  the  oral  airways. 

(2)  JP-8  fuel  droplet  evaporation  greatly  affects  deposition  concentrations  in  human  airways. 
Droplet  deposition  fractions  due  to  evaporation  decrease  with  elevated  ambient  temperatures  and 
lower  inspiratory  flow  rates. 

(3)  Assuming  idealized  inlet  velocity  profdes  and  droplet  distributions  may  greatly 
overpredict  droplet  deposition  efficiencies  for,  say,  dp>3|im,  in  the  upper  bronchial  tree. 
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(4)  Although  the  local  vapor  concentrations  may  be  affected  by  non-isothermal  flow,  the  total 
and  segmental  deposition  fractions  of  (JP-8)  fuel  vapor  in  the  upper  airways  are  basically  not 
influenced  by  variations  in  ambient  temperature,  i.e.,  263K<Tjn<310K,  where  the  variation  of 
deposition  fractions  is  within  15%. 

(5)  The  flow  rate  has  a  significant  effect  on  fuel  vapor  deposition,  i.e.,  the  higher  the  flow  rate 
the  lower  is  the  deposition  fraction;  for  example,  the  deposition  fraction  in  the  present  system, 
i.e.,  mouth  to  generation  G3,  changes  from  30%  to  60%  when  switching  from  high-level  to  low- 
level  inhalation. 

(6)  The  airway  wall  absorption  of  different  JP-8  fuel  components  may  greatly  influence  fuel 
vapor  deposition 

(7)  As  with  micro-size  droplets,  deposition  of  vapor  occurs  to  a  greater  extent  around  the 
carinal  ridges  when  compared  to  the  straight  segments  in  the  bronchial  airways;  however,  vapor 
deposition  distributions  are  much  more  uniform  along  the  airway  branches.  Such  a  near-uniform 
distribution  of  deposited  vapor  also  implies  that  more  sites  may  interact  with  cell  membranes, 
hence  there  is  a  higher  probability  of  transport  of  toxic  substances  into  tissue/blood. 

4.3  Transitions 

•  U.S.  EPA,  Health  Effects  Research  Laboratory  (Dr.  C.  S.  Kim) 

Result:  Computational  aerosol  inhalation  studies  of  focal  fuel  droplet  as  well  as 
ethanol  and  MTBE  vapor  depositions  in  a  human  upper  airway  models 
Application:  Validated  toxic  material  deposition  results  are  being  used  by  the 
EPA  for  dosimetry-and-health-effect  investigations,  involving  experimentalists, 
toxicologists,  and  pathway  modelers.  In  summary,  the  AFOSR-sponsored 
technology  can  play  a  major  role  in  analyzing  cause-and-effect  scenarios  and 
ultimately  justifying  revised  and  new  air  pollution  standards. 

•  Lovelace  Respiratory  Research  Institute  (LRRI),  Inhalation  Drug  Delivery  Center  (Dr. 
Y.S.  Cheng) 

-  Result:  Data  exchange  and  supply  of  human  cast  models 

•  UNC-CH,  Division  of  Pulmonary  and  Critical  Care  Medicine  (Dr.  J.  Donohue) 

Result:  Supply  of  project  findings  to  UNC-CH 

•  CUT  Centers  for  Health  Research  (Dr.  Julie  Kimbell) 

-  Result:  Information  exchange 
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provides  experimental  data  sets  for  lung  airway  models,  reviews  our  manuscripts,  and 
advises  on  toxic  particle  characteristics 

•  Dr.  Y.  S.  Cheng,  Director,  Inhalation  Drug  Delivery  Center,  LRRI,  Albuquerque,  NM. 
Dr.  Kim  provides  experimental  data  sets  for  lung  airway  models,  reviews  our 
manuscripts,  and  advises  on  toxic  particle  characteristics 

•  Dr.  J.  Donohue,  Chief  of  Pulmonary  Division,  UNC-CH  School  of  Medicine.  Dr. 
Donohue  advises  us  on  aerosol  characteristics,  lung  diseases,  etc. 
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